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There is much interest in how quantum systems thermalize after a sudden change, because unitary evolution 
should preclude thermalization. The eigenstate thermalization hypothesis resolves this because all 
observables for quantum states in a small energy window have essentially the same value; it is violated for 
integrable systems due to the infinite number of conserved quantities. Here, we show that when a system is 
driven by a DC electric field there are five generic behaviors: (i) monotonic or (ii) oscillatory approach to an 
infinite-temperature steady state; (iii) monotonic or (iv) oscillatory approach to a nonthermal steady state; 
or (v) evolution to an oscillatory state. Examining the Hubbard model (which thermalizes under a quench) 
and the Falicov- Kimball model (which does not), we find both exhibit scenarios (i-iv), while only Hubbard 
shows scenario (v). This shows richer behavior than in interaction quenches and integrability in the absence 
of a field plays no role. 

The classical picture for how an isolated system evolves as it is driven by a DC electric field (E) is that a current 
(j(t)) develops which subsequently creates heat due to Joule heating at a rate given by (j(t))*E. As heat is 
added to the system, the occupancy in momentum space becomes more and more uniform until the current 
decreases to zero, the occupation of all states is equally likely, and the system has reached the infinite-temperature 
limit, which is a thermal state, even though the system was driven to nonequilibrium for all earlier times. Such 
behavior should occur in any interacting system (integrable or not). This then opens the question, do driven 
quantum systems evolve in a similar way as they do under an interaction quench 1 " 4 ? There are many similarities in 
these two systems to suggest that they should. A DC field being turned on modifies the Hamiltonian at an instant 
of time, and for future times, it can be described by a time-independent Hamiltonian in the scalar potential only 
gauge. If the system thermalizes to the infinite-temperature state when it is field driven, then it also evolves into a 
steady thermal state in the infinite-time limit. 

But there are differences too, the most important being that the field- driven system has a current flowing 
through it as it thermalizes, so its energy continuously evolves as opposed to the instantaneous change in an 
interaction quench, and it is much more likely to show oscillatory behavior due to the possibility of Bloch 
oscillations. For example, driving a noninteracting single-band system with a DC field creates an oscillating 
current, and the heating varies periodically in time so that the system returns to its initial equilibrium state 
after each Bloch period given by 2n/E> resulting in no net heating; in particular, it never evolves to the 
infinite-temperature state. This occurs because there is no scattering in the system that could allow it to 
thermalize. In this work, we describe what happens for the general case when the field- driven system is 
interacting. 

We consider two interacting quantum systems: the Hubbard and the Falicov- Kimball (FK) models, which have 
been respectively shown to thermalize 5 and not to thermalize 6 under an interaction quench (the FK model has an 
infinite number of conserved quantities, but is not integrable via the Bethe ansatz in one dimension). They are 
initially in equilibrium and we study their long time behavior after a constant DC electric field is turned on at t = 
0. We then track the real-time transient behavior as the systems evolve toward a steady state and find that the 
formation of a nonequilibrium steady state density of states (DOS) is only constrained by a "causality" timescale 
(set by its Fourier transform in the time domain) that plays no further role in the relaxation of the system. At half- 
filling, an identity relates the DOS to the imaginary part of the lesser Green's function (defined and proven below), 
so in this case, it is the real part of the lesser Green's function, Re[G < (t ave ,t re i)] with t ave = (t + t')l2 the average 
time and t re i = t — t' the relative time, that determines the thermalization. It vanishes for an infinite temperature 
thermal state and is nonzero for nonequilibrium steady or oscillatory states. We use the evolution of this Green's 
function (and the current, the kinetic and the potential energies) as a function of time to describe the different 
scenarios seen for the long-time evolution of the system. 
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The equilibrium Hamiltonian describing the Hubbard 7 and FK 8 
models is 

H= - Y.ija 1 ^^- Y.ic^ c i c "+ u Y.i c hA c ^ w 

where c] 0 (c io ) are the creation (annihilation) operators for an elec- 
tron at site i with spin a, \i a is the chemical potential for the corres- 
ponding electron (which is independent of spin for the Hubbard 
model), and U is the interaction energy. The hopping integral is 
nonzero only between nearest neighbors and serves as our energy 

unit. We take J^ p = =J j ^2\fd^ for the Hubbard model, where 

d is the spatial dimension of the lattice; for the FK model, the up spin 
hopping is the same, but the down spin hopping vanishes since those 
electrons are localized and do not hop. The hopping / is used as the 
energy unit. 

The nonequilibrium case has a spatially uniform, but time varying, 
electric field, that is described by a vector-potential-only gauge, 
where E(t) = —dA(t)/dt and we have set c = e = H = 1. The hopping 

integral becomes time dependent, acquiring a phase, J^if)^]^ 
exp{ — i [Ri — Rj] -A(f) } with R { the position vector for site i on the 
lattice 9 . We start the system in equilibrium at an initial temperature, 
and then instantly turn on a constant electric field, whose spatial 
component is E, at time t = 0. The field usually points in the diagonal 
direction E = (£,£,£,. . .), or in the axial direction E = (£,0,0,. . .). The 
FK model is solved exactly on the infinite- dimensional hypercubic 
lattice, using nonequilibrium dynamical mean-field theory 1011 , while 
the Hubbard model is solved approximately using a strong- coupling 
perturbation theory in the hopping on a three-dimensional simple 
cubic lattice 12 . 

Results 

In Fig. 1, we plot the effective temperature, determined by equating 
the transient energy of the nonequilibrium system at time t to the 
equilibrium energy at temperature T in order to obtain the effective 
temperature T(t). (Note that extraction of an effective temperature by 
equating the instantaneous nonequilibrium energy to the equilib- 
rium energy at a given temperature does not necessarily imply that 
the system is in equilibrium at that instant, instead it is a convenient 
way to keep track of what an effective temperature for the system 
should be, and as the temperature approaches infinity, the system 
approaches closer and closer to a true equilibrium state.) As the time 



approaches infinity, there is a clear evolution toward an infinite- 
temperature result for cases that thermalize and a clear evolution 
to a finite temperature (indicating a nonthermal evolution, since 
the system will not be in equilibrium) for those that do not. The 
limiting behavior can be approached in a monotonic or oscillatory 
fashion (with a dynamical phase transition in between 5 ). Note that 
while the results for the FK model, where our algorithm is numer- 
ically exact, have a clearer differentiation of the different categories 
than in the Hubbard model, where our results are approximate, we 
cannot rule out the possibility that the nonthermal Hubbard model 
states are prethermalized states that subsequently evolve to the infi- 
nite-temperature limit at very long times 51314 . But our results show 
no indication of this, and such timescales are rarely relevant to any 
kind of real experiment when the interaction is large enough. One 
should also note that in some cases the effective temperature is nega- 
tive. This arises from a transient population inversion in the system. 

We begin our analysis by describing the nonequilibrium Green's 
functions which depend on two times because the system has no 
time-translation invariance. The retarded Green's function G^(t,t f ) 
and the lesser Green's functions (t,t f ) are defined via 



Gf j (t,t')=-i6(t-t f )Tre 



-0H(f->-oo) 



{c i (t),c)(t')} + /z (2) 



and 



G?(t,t') = iTre-P H ^-^c](t')ci(t) / Z, 



(3) 



where 0(t) is the unit step function, {.,.} + denotes the anticommu- 

tator, z = Tre~P Hiyt ^~ cc " > is the initial (equilibrium) partition func- 
tion and all operators are expressed in the Heisenberg representation 
with respect to the full time-dependent Hamiltonian. The local lesser 
Green's function satisfies a simple identity given by 



G<(t,t')*=-G<(f,t), 



(4) 



which follows from the definition of the lesser Green's function in Eq. 
(3) and the invariance of the trace under a cyclic reordering of its 
terms. If we express this in terms of the Wigner time coordinates, 
then we find that ReG^ (t ave , t re \) = — ReG^ (t avei — t re i) and Im 
(t aV e, t re i) = Im G^ (t ave , — t re i), or the real part of the local lesser 
Green's function is an odd function of relative time and the imagin- 
ary part is an even function of relative time. Next, we examine the 
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Figure 1 | Effective temperature as a function of time for the half-filled Falicov- Kimball and Hubbard models driven by a field (the time axis has been 
shifted so that the field is turned on at tj = 1). The colors indicate the different scenarios detailed in the text. The parameters can be found for the 
corresponding cases in the caption to Figs. 3-7. The magenta dashed line indicates the infinite temperature limit corresponding to a thermalized system. 
While we cannot rule out the possibility of the non-thermal states thermalizing on longer time scales, we see no indication of this occurring in any of the 
data we analyzed. 
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Figure 2 | Imaginary part of the local retarded Green's function as a function of relative time for various average times, (a) the Falicov- Kimball 
model and (b) the Hubbard model. The lowest curve shows the equilibrium (green) and the nonequilibrium steady state (NESS, blue) G R (t re i); The upper 
curves show G R ( t avey t re i) for successive average times after the field is switched on superposed on the steady state result (blue). Two regions are highlighted: 
(i) one in which both fand t' have the field on (red, overlapped by the blue) and the other in which fhas the field on and t' does not (orange, mixed Green's 
function). The black dots mark the causality time between the two regions. The parameters are as follows: FK model; E = 1.0, U = 3.0, T = 0.1, Hubbard 
model; U = 4^3, T = )3,E= (4^3,4^3,4^3) . Note how the Green's function becomes negligible in size when we reach the dynamic range time on the right 
hand side of each panel, indicated approximately by the gray shaded region. Panels (c) and (d) show the initial equilibrium and the final steady state DOS 
after the relative time is Fourier transformed to frequency. 



retarded Green's function under a particle-hole transformation. 
Here we assume that the lattice is a bipartite lattice, and the hopping 
is between the different sublattices only. If the chemical potential 
satisfies fi = (7/2, then the Hamiltonian is invariant under the unitary 
particle -hole transformation given by 

cJ-K-l) 6 ^ andc / ->(-l) e(/) c|, (5) 

where e({) is 1 if i is on the A sublattice and 0 if i is on the B sublattice. 
In the Hubbard model, the particle-hole transformation involves 
both spins, while for the Falicov-Kimball model, the transformation 
involves only the conduction electrons, while the localized electrons 
must be transformed via iv f —> 1 — w t . Since this is a unitary trans- 
formation, and the Hamiltonian is unchanged by it, we immediately 
see that 

-i{ Ci (t)c](t')) = -i(cj(f)c,(0)= -i(cKt)c,(f)) (6) 

where the first equality comes from the unitary particle -hole trans- 
formation (the minus signs cancel because the two operators are at 
the same site). The Hamiltonian used in evaluating the middle 
expectation value is H, and the second equality follows from the fact 
the Hamiltonian is equal to its particle-hole transformed version at 



half filling (11 = 11), which implies that we can express the expecta- 
tion value entirely in terms of tilde operators or in terms of operators 
without the tilde. The retarded Green's function then satisfies 

G f i (t,f) = -ie(t-t>)[(c i (t)J i (f)) + (c}(f) Ci (t))] 

= -iOit-t 1 ) [(cJ(() Ci (0) + {cl(t')c(t))] 

= -ie(t-t')[(cl(t)c i (t)y + (J i (t')c i (t))] 

= 9(t-t')[-G<(t/r + G<(t,t')] 

where the intermediate steps involve applying particle-hole sym- 
metry, taking the complex conjugate and applying the definition of 
the lesser Green's function. These results immediately show us that 
Re Gg(f/) = 0 at half filling and Im G?(f ,0 = 20(f - t')\m (t,t'). 
This result, which holds at half filling, tells us that the evolution of the 
local retarded Green's function as a function of its time arguments is 
identical to the evolution of the imaginary part of the local lesser 
Green's function as a function of its time arguments. 

We now examine how rapidly the local DOS approaches its 
steady-state behavior. The transient local DOS of the system is 
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Figure 3 | The oscillatory thermalization scenario. Real part of the lesser Green's function as a function of t re \ for various average times, 
(a) Falicov- Kimball, (b) Hubbard model; A grayscale is used with lighter shades indicating later average times. The inset in both cases shows the amplitude 
of G < (t rd ) as a function of average time. Panels (c) and (d) show for the FK and the Hubbard models respectively, the total energy (TE, green), the 
potential energy (PE, blue), kinetic energy (KE, red) and the current (black) as a function of time with the same parameters as in (a) and (b), respectively. 
The parameters are FK, E = 0.5, U = 0.5, T = 0.1. Hubbard, U = 4^3, E = (l,l,l)x4,/3, T = ^3. 




Figure 4 | The monotonic thermalization scenario. Real part of the lesser Green's function as a function of t rd for various average times, 
(a) Falicov- Kimball, (b) Hubbard model; A grayscale is used with lighter shades indicating later average times. The inset in both cases shows the amplitude 
of G < (t re i) as a function of average time. Panels (c) and (d) show for the FK and the Hubbard models respectively, the total energy (TE, green), the 
potential energy (PE, blue), kinetic energy (KE, red) and the current (black) as a function of time with the same parameters as in (a) and (b), respectively. 
The parameters are FK, E = 0.5, U = 1.5, T = 0.1. Hubbard, U = 4^3, E = (4^3,0,0), T = ^3. 



SCIENTIFIC REPORTS | 4:4699 | DOI: 1 0.1 038/srep04699 



4 




Figure 5 | The oscillatory nonthermal scenario. Real part of the lesser Green's function as a function of t rd for various average times, (a) Falicov- 
Kimball, (b) Hubbard model; A grayscale is used with lighter shades indicating later average times. The inset in both cases shows the amplitude of G < ( t re i) 
as a function of average time. Panels (c) and (d) show for the FK and the Hubbard models respectively, the total energy (TE, green), the potential energy 
(PE, blue), kinetic energy (KE, red) and the current (black) as a function of time with the same parameters as in (a) and (b), respectively. The parameters 
are FK, E = 2.0, U = 1.0, T = 0.1. Hubbard, U = 4^3, E = (l,l,l)x5,/3, T = ^3. 
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Figure 6 | The monotonic nonthermal scenario. Real part of the lesser Green's function as a function of t re \ for various average times, (a) Falicov- 
Kimball, (b) Hubbard model; A grayscale is used with lighter shades indicating later average times. The inset in both cases shows the amplitude of G < (t re i) 
as a function of average time. Panels (c) and (d) show for the FK and the Hubbard models respectively, the total energy (TE, green), the potential energy 
(PE, blue), kinetic energy (KE, red) and the current (black) as a function of time with the same parameters as in (a) and (b), respectively The parameters 
are FK, E = 2.0, U = 3.0, T = 0.1. Hubbard, U = 4^3, E = (5^3,0,0), T = ^3. 
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Figure 7 | The oscillatory nondecaying scenario. The lesser Green's function as a function of t re \ for various average times in (a). A grayscale is used with lighter 
shades indicating later average times. The inset shows the amplitude of G < (t rd ) as a function of average time, (b) shows the total energy (TE, green), the potential 
energy (PE, blue), kinetic energy (KE, red) and the current (black) as a function of time. The parameters are Hubbard, U = 4^3, E = (2,2,2)x4y3, T = ^3. 



determined by the local retarded Green's function, after Fourier 
transformation with respect to the relative time: A(t avei w) = — Im 
{J dt re iG R (t ave ,t re i)exip(i(Dt re i) /n}. The DOS in the frequency 
domain typically has a finite bandwidth and no singularities, and 
hence it is described by a Green function in time that has a finite 
extent in t rd that we denote as the dynamic range time t dyn . This is 
determined primarily by the inverse of the bandwidth of the density 
of states (as a function of frequency), as long as the density of states 
does not have any sharp features or structures that can give rise to 



long tails in the time domain. The fact that the dynamic range time is 
finite is important because for a fixed t ave , there always is some range 
of t rd such that t' < 0 is before the field is turned on and t > 0 is after, 
with the causality line determined by t rd = 2t ave . When both times 
have the field on, the DOS is described by the nonequilibrium steady 
state DOS. When t > 0 is after and t' < 0 before, the DOS is described 
by the mixed Green's function, interpolating between the nonequili- 
brium and original equilibrium DOS. Once the causality line passes 
the dynamic-range time (2t ave = t dyn ), the transient local DOS is 




F-K 



Figure 8 | Imaginary part of the lesser Green's function as a function of frequency for different average times in the Falicov- Kimball model at half 
filling. The system starts in the equilibrium result given by the solid green curve before the field is turned on. As the average time increases, the system 
converges to the infinite time and infinite temperature limit given by the blue dashed curve. But even before reaching that limit point, the approximation 
of the lesser Green's function by an appropriate fluctuation- dissipation theorem given by a quasi-equilibrium temperature (red curve, illustrating the 
longest simulated time) agrees very well with the corresponding transient Green's function extracted directly from the exact solution of the model (series 
of gray/black curves terminating at the red curve at the longest simulated time). The parameters are T= 0.1, U = 1.5, and E = 0.5, corresponding to a 
monotonic, thermalized case. 
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essentially given by the steady-state result because the part of the 
Green's function corresponding to relative times longer than the 
dynamic range time give essentially no contribution to the Fourier 
transform. The dynamic range time is often a short time for inter- 
acting systems, but is infinitely long for the noninteracting single 
band case, since the steady state DOS in the presence of the field is 
given by the series of delta functions describing the Wannier- Stark 
ladder, which has an infinite dynamic range time. We illustrate the 
generic situation in Fig. 2, which shows these scenarios for the 
Falicov- Kimball and Hubbard models; for both models, the system 
is in equilibrium and the field is turned on at t = 0. Hence, the 
imaginary part of the local lesser Green's function approaches it's 
steady state rapidly, just like the retarded Green's function does, and 
the dynamic evolution of the system is primarily encoded in how the 
real part of the local lesser Green's function behaves as a function of 
time (it vanishes in the infinite-temperature limit, but is nonzero in 
other situations). 

In Figs. 3-6, we show examples of the four scenarios that lead to 
thermalized or nonequilibrium steady states in the long-time limit. 
The approach to this limit (at long times) is either monotonic (over- 
damped, Figs. 4 and 6), or oscillatory (underdamped, Figs. 3 and 5), 
and the crossover between these two regimes has been called a none- 
quilibrium dynamical phase transition 5 (we primarily use the long- 
time behavior of the total energy to determine the classification). 
Each scenario is illustrated with four panels. The top panels in each 
figure show Re G < for the Falicov-Kimball model (left) and the 
Hubbard model (right), while the bottom panels in each figure show 
the current, kinetic, potential, and total energies (left) for the FK 
model and (right) for the Hubbard model. Because we cannot evolve 
these systems out to infinite time, we cannot distinguish between a 
nonthermal steady state and a transient prethermalized state that 
subsequently evolves to the infinite-temperature thermalized state 
if the prethermalized state is metastable for a long period of time. 
However, we have no evidence for that behavior either (even though 
it is often assumed to occur in the literature). Once the kinetic energy 
and current are suppressed to zero, we do not expect any further 
evolution of the system in time; the nonthermal states arise when the 
current either vanishes before the infinite-temperature state is 
reached or it oscillates equally up and down to generate no net heat- 
ing. Note that the thermalized infinite-temperature cases do not have 
simple equilibrium analogues, because they have a sharply modified 
DOS due to the driving by the field 1215 " 17 , so that even though they 
are described by thermal distribution functions, they cannot be 
described by equilibrium models without a field. 

In Fig. 7, we show the final scenario, evolution to a long-term 
oscillatory state that does not thermalize or become a steady state. 
This always occurs for a single-band noninteracting system, and here 
we found an example of this for the Hubbard model. We did not see 
such behavior in the FK model, but it might be a strong- coupling 
phenomenon, and the algorithm used to exactly solve the FK model 
cannot be accurately extended to the regime where the interaction is 
very large due to the need for too small a discretization size along the 
contour. 

Finally, we show an example of how a system that has a monotonic 
exponential approach to the infinite-temperature limit can often be 
described by a fluctuation- dissipation theorem-like picture in Fig. 8. 
This is done by first extracting an effective temperature by equating 
the instantaneous energy with the thermodynamic energy evaluated 
at equilibrium as a function of the temperature as previously shown 
in Fig. 1. Next, one takes the nonequilibrium steady state DOS (which 
can be found from short time transient calculations due to its rapid 
thermalization, as explained above, or can be solved for exactly using 
a Floquet-like theory for the Falicov-Kimball model 1617 ). Then one 
simply forms the approximate quasi- equilibrium lesser Green's func- 
tion by taking the product of the nonequilibrium steady state DOS 
with the appropriate Fermi factor for each given average time. We 



can compare this to the exact transient lesser Green's function, found 
by Fourier transforming the relevant lesser Green's function for fixed 
average time as a function of the relative time. These results are 
shown in Figure 8. The agreement between these two different 
approaches points to a quasi-equilibrium regime when the current 
is almost zero and that can be described by an effective temperature 
obtained via energy conservation. This picture is in agreement with 
the general conjecture found in the context of a one-dimensional 
system of fermions accelerated by an external field 18 . Note that the 
peaks that form in the figure are just the Wannier- Stark peaks which 
appear at multiples of the field amplitude, which are broadened by 
the interactions. 

Discussion 

We have illustrated that the thermalization problem for a field driven 
nonequilibrium system is complex, showing either monotonic or 
oscillatory approach to the thermal state or to a nonequilibrium 
steady state or evolution to an oscillating, periodic state that does 
not thermalize (with a dynamic phase transition or crossover sepa- 
rating the different regions). How the system thermalizes (or not) 
under the quench of the interaction strength does not appear to 
provide any evidence for how it will evolve when driven by a field, 
as the two models examined illustrate different behavior under a 
quench, but similar behavior when driven by a field. These results 
show that the field- driven thermalization problem has much richer 
behavior than conventional quench problems, which opens up a new 
realm for analysis of nonequilibrium behavior. The classification of 
different scenarios is simplified and clarified by examining the beha- 
vior in the time domain rather than the more commonly studied 
frequency domain. It is interesting to also study field driven systems 
in other contexts, like in the response to a pulse, where the final state 
can possibly thermalize at any temperature since the final Hamil- 
tonian is the same as the initial Hamiltonian the system had in 
equilibrium, and in the long time limit no current flows 19,20 . 

Methods 

The strong-coupling perturbation theory for the Hubbard model uses a self- 
consistent second-order expansion for the self-energy, which includes an infinite 
class of diagrams, but is truncated 12 . The exact solution for the FK model uses a 
discretized version of nonequilibrium dynamical mean-field theory which was 
extrapolated to the continuum limit 10,11 . The Kadanoff-Baym-Keldysh technique is 
employed to solve for the Green's functions. There are two Green's functions to 
determine. The retarded Green's function Gfj(t,t f ) and the lesser Green's function 
G~ (t,t') defined in Eqs. (2) and (3). The corresponding momentum-dependent 
Green's functions are formed by taking the spatial Fourier transform, since the 
Green's functions depend only on the spatial difference of the position variables 
(due to translational in variance of the system). These Green's functions are cal- 
culated with the numerically exact dynamical mean-field theory approach for the 
Falicov-Kimball model and with the self-consistent second-order strong- coupling 
expansion for the Hubbard model. 

The Green's functions can be employed to calculate the kinetic energy, potential 
energy, and current in the presence of a driving field 20 . The bandstructure satisfies 

s(k = A(t)) = — i^J^j cos (^' ~Mt)) i n the limit of d approaching infinity, 

with A(t) the component of the vector potential, which is nonzero only for positive 
times. The expectation value for the kinetic energy then becomes 

£Ut) = J2k <k-A(t))(cl(t)c k (t)) = ->J2 k <k-A(t))Gt (t,t) (8) 

which is determined from the momentum-dependent lesser Green's function (and is 
per spin, so must be multiplied by two for the Hubbard model). Similarly, if we define 

the band velocity via v( k — A (t)) = — {^J^j Si-i sm (^' — ^M) also in the limit of 
d going to infinity, then the current becomes 

j(t) = J2 k v(k-A(t))(cl(t)c k (t)) = -iJ2 k v(*"A(f))G fc < (M) (9) 

which is also per spin and must be multiplied by two for the Hubbard model. The 
potential energy is proportional to the average double occupancy and requires the 
time derivative of the Green's function to evaluate it. It satisfies 
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17, 

dG£ (t, t') , 

: df 



V(0 = ^I^(«<r(*HW) 

3(7 



L7 



(10) 



" \t^t + — (rt T (t) + »i W) - - (0- 

The total energy is just the sum of the kinetic plus potential energies. 
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